* Independently-generated stationary outcome and unitroot predictor DGP

// Basics
clear all
set matsize 11000
parallel setclusters 8

capture program drop NG
program define NG, rclass

// Monte Carlo Basics
version 12.1
syntax [, obs(integer 1) phi_y(real 0.5) phi_x(real 0.5) ]
drop _all
set obs `obs'

// Time Series Basics
egen time = fill(0 1 2)
tsset time

// Data Generating Process

* Generate predictor
gen x = rnormal()
replace x = `phi_x' * L.x + rnormal() if time > 0

* Generate outcome
gen y = rnormal()
replace y = `phi_y' * L.y + rnormal() if time > 0

// Model Estimation

* Estimate static model
reg y x
return scalar b1_stat = _b[x]
scalar b1_stat_test = (_b[x] - 0) / _se[x]
return scalar b1_stat_t1 = 2*ttail(e(df_r),abs(b1_stat_test)) < 0.05

* Estimate Differences model 
reg D.y D.x
return scalar b1_diff = _b[D1.x]
scalar b1_diff_test =  (_b[D1.x] - 0) / _se[D1.x]
return scalar b1_diff_t1 = 2*ttail(e(df_r),abs(b1_diff_test)) < 0.05

* Estimate partial adjustment 
reg y L.y x
return scalar phi_pa = _b[L1.y]
scalar phi_pa_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi_pa_t1 = 2*ttail(e(df_r),abs(phi_pa_test)) < 0.05
return scalar b1_pa = _b[x]
scalar b1_pa_test = (_b[x] - 0) / _se[x]
return scalar b1_pa_t1 = 2*ttail(e(df_r),abs(b1_pa_test)) < 0.05
return scalar lrm_pa = _b[x]/(1-_b[L1.y])
nlcom (1e+3)* (_b[x]/(1-_b[L1.y]))
matrix define var_lrm_pa = r(V)
scalar var_lrm_pa_se  = sqrt(var_lrm_pa[1,1])/(1e+3)
scalar lrm_pa_test = ((_b[x]/(1-_b[L1.y])) - 0) /var_lrm_pa_se
return scalar lrm_pa_t1 = 2*ttail(e(df_r),abs(lrm_pa_test)) < 0.05

scalar phi_pa_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_pa_power = 2*ttail(e(df_r),abs(phi_pa_test_power)) < 0.05

* Estimate dead start 
reg y L.y L.x
return scalar phi_ds = _b[L1.y]
scalar phi_ds_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi_ds_t1 = 2*ttail(e(df_r),abs(phi_ds_test)) < 0.05
return scalar b2_ds = _b[L1.x]
scalar b2_ds_test = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_ds_t1 = 2*ttail(e(df_r),abs(b2_ds_test)) < 0.05

scalar phi_ds_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_ds_power = 2*ttail(e(df_r),abs(phi_ds_test_power)) < 0.05

* Estimate GECM 
reg D.y L.y D.x L.x
return scalar phi_gecm = _b[L1.y]
local gecm_true_phi = `phi_y' - 1
scalar phi_gecm_test = (_b[L1.y] - `gecm_true_phi') / _se[L1.y]
return scalar phi_gecm_t1 = 2*ttail(e(df_r),abs(phi_gecm_test)) < 0.05
return scalar b1_gecm  = _b[D1.x]
scalar b1_gecm_test = (_b[D1.x] - 0) / _se[D1.x]
return scalar b1_gecm_t1 = 2*ttail(e(df_r),abs(b1_gecm_test)) < 0.05
return scalar b2_gecm  = _b[L1.x]
scalar b2_gecm_test = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_gecm_t1 = 2*ttail(e(df_r),abs(b2_gecm_test)) < 0.05
return scalar lrm_gecm = -_b[L1.x] / _b[L1.y]
nlcom (-_b[L1.x] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm = r(V)
scalar var_lrm_gecm_se  = sqrt(var_lrm_gecm[1,1]) / (1e+3)
scalar lrm_gecm_test = ((-_b[L1.x] / _b[L1.y]) - 0) /var_lrm_gecm_se
return scalar lrm_gecm_t1 = 2*ttail(e(df_r),abs(lrm_gecm_test)) < 0.05

scalar phi_gecm_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_gecm_power = 2*ttail(e(df_r),abs(phi_gecm_test_power)) < 0.05

* Estimate ADL 
reg y L.y x L.x
return scalar phi_adl = _b[L1.y]
scalar phi_adl_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi_adl_t1 = 2*ttail(e(df_r),abs(phi_adl_test)) < 0.05
return scalar b1_adl = _b[x]
scalar b1_adl_test = (_b[x] - 0) / _se[x]
return scalar b1_adl_t1 = 2*ttail(e(df_r),abs(b1_adl_test)) < 0.05
return scalar b2_adl = _b[L1.x]
scalar b2_adl_test = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl_t1 = 2*ttail(e(df_r),abs(b2_adl_test)) < 0.05
return scalar lrm_adl = (_b[x] + _b[L1.x]) /(1-_b[L1.y])
nlcom ((_b[x] + _b[L1.x]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl = r(V)
scalar var_lrm_adl_se  = sqrt(var_lrm_adl[1,1]) / (1e+3)
scalar lrm_adl_test = (((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - 0) /var_lrm_adl_se
return scalar lrm_adl_t1 = 2*ttail(e(df_r),abs(lrm_adl_test)) < 0.05

scalar phi_adl_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_adl_power = 2*ttail(e(df_r),abs(phi_adl_test_power)) < 0.05

* Estimate ADL (2,2) 
reg y L.y L2.y x L.x L2.x
return scalar phi1_adl22 = _b[L1.y]
scalar phi1_adl22_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi1_adl22_t1 = 2*ttail(e(df_r),abs(phi1_adl22_test)) < 0.05
return scalar phi2_adl22 = _b[L2.y]
scalar phi2_adl22_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl22_t1 = 2*ttail(e(df_r),abs(phi2_adl22_test)) < 0.05
return scalar b1_adl22 = _b[x]
scalar b1_adl22_test = (_b[x] - 0) / _se[x]
return scalar b1_adl22_t1 = 2*ttail(e(df_r),abs(b1_adl22_test)) < 0.05
return scalar b2_adl22 = _b[L1.x]
scalar b2_adl22_test = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl22_t1 = 2*ttail(e(df_r),abs(b2_adl22_test)) < 0.05
return scalar b3_adl22 = _b[L2.x]
scalar b3_adl22_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl22_t1 = 2*ttail(e(df_r),abs(b3_adl22_test)) < 0.05
return scalar lrm_adl22 = (_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y]))*(1e+3)
matrix define var_lrm_adl22 = r(V)
scalar var_lrm_adl22_se  = sqrt(var_lrm_adl22[1,1]) / (1e+3)
scalar lrm_adl22_test = (((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_se
return scalar lrm_adl22_t1 = 2*ttail(e(df_r),abs(lrm_adl22_test)) < 0.05


scalar phi1_adl22_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl22_power = 2*ttail(e(df_r),abs(phi1_adl22_test_power)) < 0.05

* Estimate ADL (3,3) 
reg y L.y L2.y L3.y x L.x L2.x L3.x
return scalar phi1_adl33 = _b[L1.y]
scalar phi1_adl33_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi1_adl33_t1 = 2*ttail(e(df_r),abs(phi1_adl33_test)) < 0.05
return scalar phi2_adl33 = _b[L2.y]
scalar phi2_adl33_test = (_b[L2.y] - 0) / _se[L2.y]
return scalar phi2_adl33_t1 = 2*ttail(e(df_r),abs(phi2_adl33_test)) < 0.05
return scalar phi3_adl33 = _b[L3.y]
scalar phi3_adl33_test = (_b[L3.y] - 0) / _se[L3.y]
return scalar phi3_adl33_t1 = 2*ttail(e(df_r),abs(phi3_adl33_test)) < 0.05
return scalar b1_adl33 = _b[x]
scalar b1_adl33_test = (_b[x] - 0) / _se[x]
return scalar b1_adl33_t1 = 2*ttail(e(df_r),abs(b1_adl33_test)) < 0.05
return scalar b2_adl33 = _b[L1.x]
scalar b2_adl33_test = (_b[L1.x] - 0) / _se[L1.x]
return scalar b2_adl33_t1 = 2*ttail(e(df_r),abs(b2_adl33_test)) < 0.05 
return scalar b3_adl33 = _b[L2.x]
scalar b3_adl33_test = (_b[L2.x] - 0) / _se[L2.x]
return scalar b3_adl33_t1 = 2*ttail(e(df_r),abs(b3_adl33_test)) < 0.05
return scalar b4_adl33 = _b[L3.x]
scalar b4_adl33_test = (_b[L3.x] - 0) / _se[L3.x]
return scalar b4_adl33_t1 = 2*ttail(e(df_r),abs(b4_adl33_test)) < 0.05
return scalar lrm_adl33 = (_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) * (1e+3)
matrix define var_lrm_adl33 = r(V)
scalar var_lrm_adl33_se  = sqrt(var_lrm_adl33[1,1]) / (1e+3)
scalar lrm_adl33_test = (((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) - 0) /var_lrm_adl33_se
return scalar lrm_adl33_t1 = 2*ttail(e(df_r),abs(lrm_adl33_test)) < 0.05

scalar phi1_adl33_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl33_power = 2*ttail(e(df_r),abs(phi1_adl33_test_power)) < 0.05

* F-test, ADL(2,2)
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = `phi_y'
return scalar ftest_y_adl22_t1 = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test L.y + L2.y = 0
return scalar ftest_y_adl22_power = r(p) < 0.05
reg y L.y L2.y x L.x L2.x
test x + L.x + L2.x = 0
return scalar ftest_x_adl22_t1 = r(p) < 0.05

* F-test, ADL(3,3)
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = `phi_y'
return scalar ftest_y_adl33_t1 = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test L.y + L2.y + L3.y = 0
return scalar ftest_y_adl33_power = r(p) < 0.05
reg y L.y L2.y L3.y x L.x L2.x L3.x
test x + L.x + L2.x + L3.x= 0
return scalar ftest_x_adl33_t1 = r(p) < 0.05


end
